


rm(list=ls())
set.seed(1)

library(mediation)
library(sandwich)
library(Zelig)
library(VGAM)

setwd("C:/Users/dtingley/Documents/My Dropbox/Tingley AnxietyStudies/BargainingProtocolData/Data/JournalReplication/Distribution")
setwd("C:/Users/dtingley/Dropbox/Tingley AnxietyStudies/BargainingProtocolData/Data/JournalReplication/Distribution")

load("IOReplicationDataPhysio.RData")
source("ci-plot-new.R")

setwd("../../../Paper/JournalFiles")

###
#Descriptive Stats
###

pdf("OfferDistributionRound0.pdf", onefile=FALSE, paper="special")
par(mfrow=c(2,1),mar=c(4,3.75,2,1), oma=c(.5,.1,.2,.2), cex.axis=.8, cex.main=.8, cex.lab=.8)
hist(anxiety.0$AOffer[anxiety.0$roundId==0], freq=FALSE, main="Small Shift in Power",xlab="Offer")
abline(v=mean(anxiety.0$AOffer[anxiety.0$roundId==0]),lty=2)
hist(anxiety.1$AOffer[anxiety.1$roundId==0],freq=FALSE,  main="Big Shift in Power",xlab="Offer")
abline(v=mean(anxiety.1$AOffer[anxiety.1$roundId==0]),lty=2)
dev.off()



####
#Zelig plots
####



#Only compare large and big shift with offer changes
#estimate by interacting treatment w offer.
ProbReject<-zelig(BDecision ~treat.nofix*AOffer , model="probit", data=subset(anxiety,roundId==0))
BigShift<-setx(ProbReject,treat.nofix=1,AOffer=c(seq(0,10,by=.1)))
SmallShift<-setx(ProbReject,treat.nofix=0,AOffer=c(seq(0,10,by=.1)) )
s.out1<-sim(ProbReject,x=SmallShift, x1=BigShift, num = 10000)

pdf("RejectPlot.pdf", onefile=FALSE, paper="special")
ci.plot(s.out1,var="AOffer", ci=c(95), leg=0, xlab="Offer", ylab="Probability of Rejection" )
text(x=7.5,y=0.4,label="Big Shift in Power")
text(x=3.5,y=0.4,label="Small Shift in Power")
dev.off()


skin.2<-zelig(SkSDRecOffer1B_mean_r  ~ AOffer , model="ls", robust=TRUE, data=subset(anxiety.2,roundId==0))
skin<-zelig(SkSDRecOffer1B_mean_r  ~ AOffer*treat.nofix , model="ls", robust=TRUE, data=subset(anxiety,roundId==0))
BigShift<-setx(skin,treat.nofix=1,AOffer=c(seq(0,10,by=.1)) )
SmallShift<-setx(skin,treat.nofix=0,AOffer=c(seq(0,10,by=.1)))
s.out1<-sim(skin,x=SmallShift, x1=BigShift , num = 10000)

pdf("SkinPlot.pdf", onefile=FALSE, paper="special")
ci.plot(s.out1,var="AOffer", ci=c(90), leg=0, xlab="Offer", ylab="Average Skin Conductance Score" )
text(x=8,y=0.6,label="Small Shift in Power")
text(x=3,y=0.75,label="Big Shift in Power")
dev.off()



NoOfferShift<-setx(skin.2,AOffer=c(seq(0,10,by=.1)))
s.out1<-sim(skin.2,x=NoOfferShift, num = 10000)
summary(s.out1)
pdf("SkinPlotNoChange.pdf", onefile=FALSE, paper="special")
ci.plot(s.out1,var="AOffer", ci=c(90), leg=0, xlab="Offer", ylab="Average Skin Conductance Score" )
text(x=8,y=0.6,label="No Offer Change Allowed")
dev.off()





###
#Mediation
###



## Mediation model
med.fit.2 <- lm(SkSDRecOffer1B_mean_r ~ AOffer, data=subset(anxiety.2,roundId==0))
med.fit.1 <- lm(SkSDRecOffer1B_mean_r ~ AOffer, data=subset(anxiety.1,roundId==0))
med.fit.0 <- lm(SkSDRecOffer1B_mean_r ~ AOffer, data=subset(anxiety.0,roundId==0))

summary(med.fit.0)
summary(med.fit.1)
summary(med.fit.2)

## Outcome model. With   interaction
out.fit.2 <- glm(BDecision ~ SkSDRecOffer1B_mean_r*AOffer, data=subset(anxiety.2,roundId==0), family = binomial("probit"))
out.fit.1 <- glm(BDecision ~ SkSDRecOffer1B_mean_r*AOffer, data=subset(anxiety.1,roundId==0), family = binomial("probit"))
out.fit.0 <- glm(BDecision ~ SkSDRecOffer1B_mean_r*AOffer, data=subset(anxiety.0,roundId==0), family = binomial("probit"))

summary(out.fit.0)
summary(out.fit.1)
summary(out.fit.2)

## Computing the ACME etc.
l<-5
h<-10
sims<-1000
med.out.2 <- mediate(med.fit.2, out.fit.2, treat = "AOffer", mediator = "SkSDRecOffer1B_mean_r",cluster = anxiety.2$ID, sims=sims, long=FALSE, dropobs=TRUE, control.value=l, treat.value=h)
med.out.1 <- mediate(med.fit.1, out.fit.1, treat = "AOffer", mediator = "SkSDRecOffer1B_mean_r",cluster = anxiety.1$ID, sims=sims, long=FALSE, dropobs=TRUE, control.value=l, treat.value=h)
med.out.0 <- mediate(med.fit.0, out.fit.0, treat = "AOffer", mediator = "SkSDRecOffer1B_mean_r",cluster = anxiety.0$ID, sims=sims, long=FALSE, dropobs=TRUE, control.value=l, treat.value=h)

summary(med.out.2)
summary(med.out.1)
summary(med.out.0)

l<--.6
h<-.1
pdf("MediateShift510.pdf", onefile=FALSE, paper="special", height=4)
par(mfrow=c(2,1),mar=c(1.5,2,2,1), oma=c(3,3,.2,.2), cex.axis=.8, cex.main=.8, cex.lab=.8)
plot(med.out.1 , main="Big Shift in Power",xlim=c(l,h))
plot(med.out.0, main="Small Shift in Power",xlim=c(l,h))
title(xlab="Effects", outer=TRUE, line=1)
dev.off()


pdf("MediateShift510BigShiftNoOfferShift.pdf", onefile=FALSE, paper="special", heigh=4)
par(mfrow=c(1,1),mar=c(1.5,2,2,1), oma=c(3,3,.2,.2), cex.axis=.8, cex.main=.8, cex.lab=.8)
plot(med.out.2 , main="Big Shift in Power, No Change in Offer Allowed")
title(xlab="Effects", outer=TRUE, line=1)
dev.off()


























